#ifndef __SURFACE__
#define __SURFACE__

#include "primitive.h"
#include "dispatcher.h"

struct DerivativesDim2 : public Derivatives {
	DerivativesDim2(int d); 
	boost::multi_array<Vec3n, 2> D;
	boost::multi_array<Vec4n, 2> Dw;
	int d;
};

class Surface : public Primitive {
public:
	Surface() {}
	Surface(const std::vector< std::vector<Vec4n> >& Pw, int p, int q, const Knotvector& U, const Knotvector& V);
	Surface(const Array2V4n& Pw, int p, int q, const Knotvector& U, const Knotvector& V);
	Surface(const Surface& surface);

	static Surface WeightedLeastSquares(const std::vector< std::pair<Vec2n, Vec3n> >& p, size_t degree, const std::vector<nfloat>& thetad, bool usePowerBasis = false);
	static Surface LeastSquares(const std::vector< std::pair<Vec2n, Vec3n> >& p, size_t degree, bool usePowerBasis = false);

	virtual void Derivative(nfloat* param, Derivatives* pd) const;
	Vec3n Evaluate(const boost::array<nfloat, 2>& param) const;
	Vec3n Evaluate(nfloat u, nfloat v) const;
	virtual void Refine();
	void RefineU();
	void RefineV();
	virtual void Refine(const std::vector<Knotvector>& X);
	Matrix Jacobian(const boost::array<nfloat, 2>& param) const;

	virtual void Save(const std::string& filename, const std::string& attribute);
	void Save(const std::string& filename);
	
	Knotvector& GetKnotvectorU() {return U;}
	Knotvector& GetKnotvectorV() {return V;}
	const Knotvector& GetKnotvectorU() const {return U;}
	const Knotvector& GetKnotvectorV() const {return V;}
	size_t GetDegreeU() const {return p;}
	size_t GetDegreeV() const {return q;}

	void DegreeElevateU(int t);
	void DegreeElevateV(int t);

	std::vector<Vec4n> GetRow(int i) const;
	std::vector<Vec4n> GetColumn(int i) const;
	std::vector<Vec4n> SetRow(int i);
	std::vector<Vec4n> SetColumn(int i);
	void SetRow(int i, const std::vector<Vec4n>& pw);
	void SetColumn(int i, const std::vector<Vec4n>& pw);

	void Curvature(const boost::array<nfloat, 2>& param, nfloat& k1, nfloat& k2, Vec3n& d1, Vec3n& d2) const;
	void Curvature(const boost::array<nfloat, 2>& param, nfloat& k1, nfloat& k2, Vec2n& d1, Vec2n& d2) const;
	void Curvature(const boost::array<nfloat, 2>& param, boost::array<nfloat, 6>& C, nfloat& k1, nfloat& k2) const;

	void CurvatureFirstDerivative(const boost::array<nfloat, 2>& param, const bool& is_k1, nfloat& ku, nfloat& kv) const;
	void CurvatureSecondDerivative(const boost::array<nfloat, 2>& param, const bool& is_k1, nfloat& kuu, nfloat& kuv, nfloat& kvv) const;

	bool fD(Vec2n& c, const Vec2n& u, nfloat h, nfloat sign) const;
	std::vector<Vec2n> Trace(const Vec2n& u0, nfloat h) const;
	Vec3n Normal(const boost::array<nfloat, 2>& param) const;
	Vec3n Normal(nfloat u, nfloat v) const;
	std::pair<Vec3n, Vec3n> Partials(const boost::array<nfloat, 2>& param) const;
	std::pair<Vec3n, Vec3n> Partials(const nfloat u, const nfloat v) const;
	
	size_t Size0() {return Pw.size();}
	size_t Size1() {return Pw[0].size();}
	size_t Size0() const {return Pw.size();}
	size_t Size1() const {return Pw[0].size();}

	Vec4n& operator()(int i, int j) {
		if(i < 0) i = Pw.size() + i;
		if(j < 0) j = Pw[0].size() + j;
		return Pw[i][j];
	}

	const Vec4n& operator()(int i, int j) const {
		if(i < 0) i = Pw.size() + i;
		if(j < 0) j = Pw[0].size() + j;
		return Pw[i][j];
	}

	void Transpose();
	void FlipU();
	void FlipV();

	virtual void Dispatch(Dispatcher* d) {
		d->Visit(this);
	}
protected:
	void SurfaceDerivsAlg(nfloat* param, DerivativesDim2* pd) const;
	size_t p;
	size_t q;
	size_t n;
	size_t m;
	Knotvector U;
	Knotvector V;
	std::vector< std::vector<Vec4n> > Pw;
	
	void RefineKnotVectSurface(int n, int p, const Knotvector& U, 
				   int m, int q, const Knotvector& V,
				   const std::vector< std::vector<Vec4n> >& Pw, const Knotvector& X, int r, int dir, 
				   Knotvector& Ubar, Knotvector& Vbar, std::vector< std::vector<Vec4n> >& Qw);
	void binomialCoef(Matrix& Bin) const;
};

#endif //__SURFACE__
